A boundary migration model for imaging within volumetric scattering media

Effectively imaging within volumetric scattering media is of great importance and challenging especially in macroscopic applications. Recent works have demonstrated the ability to image through scattering media or within the weak volumetric scattering media using spatial distribution or temporal characteristics of the scattered field. Here, we focus on imaging Lambertian objects embedded in highly scattering media, where signal photons are dramatically attenuated during propagation and highly coupled with background photons. We address these challenges by providing a time-to-space boundary migration model (BMM) of the scattered field to convert the scattered measurements in spectral form to the scene information in the temporal domain using all of the optical signals. The experiments are conducted under two typical scattering scenarios: 2D and 3D Lambertian objects embedded in the polyethylene foam and the fog, which demonstrate the effectiveness of the proposed algorithm. It outperforms related works including time gating in terms of reconstruction precision and scattering strength. Even though the proportion of signal photons is only 0.75%, Lambertian objects located at more than 25 transport mean free paths (TMFPs), corresponding to the round-trip scattering length of more than 50 TMFPs, can be reconstructed. Also, the proposed method provides low reconstruction complexity and millisecond-scale runtime, which significantly benefits its application.


Supplementary Note 1: Equipment details and imaging scenarios
As shown in Supplementary Fig. 1a, our hardware imaging system consists of a 20fs pulsed laser (Coherent Vitara-T), a 32 × 32 array of single-photon avalanche diode detector (Photon Force PF32, 6.95µ ⌀ active area, 50µ pixel pitch), several controlling electronics and optical elements used for scanning and focusing. Supplementary Fig.  1b shows the optical path of the confocal configuration used in our system. The laser hand emits a collimated light with a wavelength of 808nm at an 80MHz repetition rate and 580mW average power. Then the light passes through a polarized beam splitter cube (Thorlabs PBS255) and is incident into the scene after being reflected by a 2-axis scanning galvanometer (Thorlabs GVS012) controlled with a data acquisition device (NI-DAQ USB-6343). Light returning from the scenario along the same optical path is focused to the detector by a lens with 8mm focal length and 1.4 maximal f-number. The delayer is used to reshape the synchronization signal output from the laser into a standard Transistor-Transistor Logic (TTL) signal, which is then used as the acquisition trigger signal of the single-photon avalanche diode (SPAD). Since the trigger signal frequency of the detector is 80MHz and its temporal resolution is 55ps, the temporal dimension of the data collected by the SPAD is 230. Therefore, the SPAD array captures the data with a size of 32 × 32 × 230 for each scanning point. We select several pixels in the array that are less affected by the initial bright light caused by the internal reflection of the confocal optical path. Then these pixels are averaged to generate the measurement of the current scanning point. By controlling the angle of the galvanometer, the system scans a grid of 64 × 64 points on the front surface of the media so that the detector collects 3D data: two spatial dimensions with the size of the scanning grid and one temporal dimension of 230 bins.
We evaluate the reconstruction performance of our algorithm in two typical scattering scenarios: 2D and 3D Lambertian objects embedded in the polyethylene foam and the fog. The two scattering imaging scenarios are shown in Supplementary Fig. 2-3. The characteristics, including scattering, absorption, and propagation angle of two scattering media, are quite different and are estimated in Supplementary Note 2. Polyethylene foam is a homogeneous static isotropic scattering medium while the fog is a dynamic forward scattering medium [1][2][3] . They produce different temporal responses as shown in Fig. 3b and 4b in the main text. The former is affected by the backscattered light so that signal photons are highly coupled with background photons. The latter is affected by the time variability of the fog and thus the temporal distribution of signal photons is superimposed with dynamic noise. Both cases exacerbate the reconstruction difficulty while we demonstrate the ability of our method in these two scattering scenarios. Supplementary Fig. 2 Photograph and layout of the polyethylene foam scattering scenario. a Confocal imaging system captures the photons reflected from the foam scenario. b The total thicknesses of the foam range from 12cm to 25cm. The object is embedded in the foam and its depth is defined as the distance from the front surface of the scattering media to the object. Supplementary Fig. 3 Photograph and layout of the fog scattering scenario. a An artificial fog chamber with a size of 50cm × 50cm × 150cm provides the dynamic scattering environment. b The object is embedded in the fog and its depth range is defined as the distance range from the front surface of the scattering media to the object. In this section, the scattering characteristics of the polyethylene foam and the fog are estimated using the solution of the diffusion equation for the slab media and Beer-Lambert Law, respectively. The former provides the optimal reduced scattering coefficient ! " and absorption coefficient # of the polyethylene foam. The reduced scattering coefficient is given by ! " = ! (1 − ), where ! is the scattering coefficient and ∈ [−1,1] is a parameter of Henyey-Greenstein phase function that describes the propagation angle of the scattering particles. Positive indicates the forward scattering and negative indicates the backward scattering. Polyethylene foam is a homogeneous static isotropic scattering medium with = 0, where the photon scatters in all angles with equal probability, while the fog is a dynamic forward scattering medium 1-3 with ≈ 0.85. The Beer-Lambert Law gives the estimation of the transport coefficient $ of the fog, which is defined as $ = ! + # .
Estimating the scattering characteristic of polyethylene foam. We estimate the reduced scattering coefficient ! " and absorption coefficient # of the polyethylene foam by optimizing the residual error between the measured temporal response and the theoretical solution of the diffusion equation 4, 5 (Eq. (1) in the main text) for the slab media which is given by: (1) where ( , ) represents the power that photons enter at the origin of the coordinate system at = 0, and exit from the far boundary of the scattering media at time per unit of area and unit of time. is the distance from the exit point to the plane = 0 as shown in Supplementary Fig. 4. Additionally, diffusion coefficient is given by = (3( ! " + # )) (& , is the speed of the light and is the thickness of the scattering media. Supplementary Equation (1) is obtained with two assumptions: (1) incident light is initially scattered isotropically at a distance of % = 1/ ! " in the scattering media 6 ; (2) the average diffuse intensity is assumed to be zero at two extrapolated flat surfaces located at a distance, ) , away from either boundary of the scattering media. To satisfy these boundary conditions, an infinite number of positive and negative dipole sources can be virtually placed in the scenario at the following positions 5 : (2) In general, 7 dipoles are sufficient to maintain the truncation error at a negligible level 5 . The extrapolated distance ) is given by: where * is the angle of incidence light, ( * ) is the reflection coefficient calculated by the Fresnel reflection coefficients 7 . For a given * , the reflection coefficient ( * ) for unpolarized light is given by: (4) where = ' / & is the relative refractive index and is defined as the ratio between the refractive index of the scattering media and that of the external media, as shown in Supplementary Fig. 4. $ = arcsin [ • sin ( * )] is the refracted angle.
Therefore, there are three parameters of scattering media to be estimated: reduced scattering coefficient ! " , absorption coefficient # and relative refractive index . The pulsed light is used to illuminate one side of the polyethylene foam, and the SPAD array is used to capture the temporal response ( ) at the other side. For simplicity, the detection point is set at = 0 then the distance is equal to 0. We measure 9 different thicknesses, , of the polyethylene foam, ranging from 2cm to 10cm with an increment of 1cm. Given a fixed value of the relative refractive index % , we optimize the ! " and # to minimize the following objective function: The Levenberg-Marquardt algorithm is used to iteratively solve this non-linear optimization in MATLAB. There are a couple of factors that need to be taken into account before the optimization. The first one is the scale of the measured temporal response which depends on the power of the incident light. So, the measurements and the theoretical values are normalized to a maximum of 1. The second one is the experimental delay between triggering the laser and the photons being captured by the detector, which depends on the optical path between the laser and the detector and the time jitter of the system. So, we need to shift the measurements to produce a good fit of the model. First, we calibrate the distance between the laser and detector. Then the offset value of each measurement is also regarded as an optimization parameter. In this case, the number of the optimization parameter in our experiment is 11, including the reduced scattering coefficient, absorption coefficient, and 9 offsets of the measurements with different thicknesses of the polyethylene foam.
Empirically given the relative refractive index across a range from 1.12 to 1.32 with an increment of 0.01, the initializations of ! " from 2.0cm (& to 3.0cm (& with an increment of 0.1cm (& and # equals 0, a scatterplot of the mean squared errors (MSE) between the model and all measurements is shown in Supplementary Fig. 5a. It can be found that the relative refractive index of 1.23 produces the minimal MSE with the initializations of ! " = 2.6cm (& and # = 0cm (& , in which case the optimal parameters are given by ! "   3.3348 × 10 (' cm (& . As shown in Supplementary Fig. 6, these parameters produce a good fit to the measured temporal responses in all thicknesses of the polyethylene foam. To quantify the uncertainty of the estimated parameters, we analyze the MSE against the reduced scattering coefficient and the absorption coefficient for all optimizations, as shown in Supplementary Fig. 5b-c. The black dashed lines indicate the optimal parameters and the red dashed lines show the boundary of the range containing 80% of all estimated values. Therefore, the reduced scattering coefficient ! " and the absorption coefficient # of the polyethylene foam used in our experiments are estimated as: ! " = 3.1377 ± 0.35cm (& and # = 3.3348 × 10 (' ± 2.6 × 10 (+ cm (& .
Supplementary Fig. 6 The fitting results of the measured temporal responses and the theoretical values. We capture 9 measurements of the polyethylene foam with different thicknesses, ranging from 2cm to 10cm with an increment of 1cm. The orange lines are the theoretical values and the blue dashed lines are the measurements. In all cases, the parameters, ! " = 3.1377cm (& and # = 3.3348 × 10 (' cm (& , produce a good fit between the measurements and the theoretical model. Estimating the scattering characteristic of the fog. The light propagates in the fog is attenuated by the scattering and absorption of the scattering particles. Both effects can be considered as a single transport coefficient $ , which is defined as the sum of the scattering coefficient ! and the absorption coefficient # . According to the Beer-Lambert Law, if the light propagates a distance with the transport coefficient $ , the optical power is given by: where % is the initial power of the incident light and $ gives the multiples of the mean free paths (MFP, defined as 1/ $ ), also known as optical thickness . In this case, the transport coefficient can be estimated by optical power measurements. The transmittance measurements are retrieved as shown in Supplementary Fig. 7, the laser emits the light with different initial optical powers into the scattering media on the one side of the fog chamber while a power meter detects the outgoing power of the light through the fog, r = 50cm, on the other side of the scattering media. We measure 10 outgoing optical powers at different initial optical powers, as shown in Supplementary Table 1 In this section, we analyze the validity range of our algorithm. The proposed BMM algorithm is based on the diffusion equation. It treats light propagating through the scattering media as a random walk 4 so that photons diffused inside the scattering media follow some defined probability rules, which can be regarded as the quantitative indications to determine whether the diffusion approximation breaks down. In particular, the average arrival time of photons is one of the quantitative indications. The average arrival time, described as the mean time required by a photon going through the scattering media, is computed by 8 where ( ′) is the temporal response measured experimentally. The theoretical prediction of the average arrival time for a photon travels a distance is proportional to ' /6 , where is the diffusion coefficient, and is the speed of light [8][9][10] . Considering the influence of the experimental factors such as the time duration of the incident pulse, the linear relationship between measured temporal displacement and theoretical average arrival time is given by 11 : where is a proportionality factor that depends on the short pulse duration and so on. We analyze the linear relationship between the measured average arrival time and the theoretical value to verify that all our experiments are Supplementary Fig. 8 The linear relationship between the measured average arrival time and the theoretical value. The relationship between the peak shift of the measured temporal response Δ and theoretical average arrival time ' /6 is fitted as a linear function. The proportionality factor is linearly fitted as 0.4812. within the validity range of the diffusion equation. We measure the temporal responses in the polyethylene foam with 9 different thicknesses ranging from 2cm to 10cm with an increment of 1cm. The measured temporal responses are shown in Supplementary Fig. 6. Then the peak shifting time of each temporal response is approximated as the measured average arrival time Δ here 11 . The fitted result indicates that measured Δ is a linear function of the theoretical value ' /6 , depicted by the orange curve in Supplementary Fig. 8, where the proportionality factor is 0.4812. It demonstrates that all the experiments are within the validity range of the diffusion equation.
Several studies show that the diffusion equation models the scattered light more accurate with the increase of the imaging depth 8,12,13 . Therefore, the validity range of our algorithm is theoretically unbounded. However, in practical experiments, the imaging quality is limited by the signal-to-background ratio (SBR) of measurement which depends on the characteristics of the physical setup, including the incident power of the laser, the sensitivity of the detector, the albedo of the object surface and so on. To quantitatively analyze the proportion of signal photons and the scattering strength in our experiments, we calculate the SBR of measurements. The SBR is defined as the ratio of signal photon intensity to background photon intensity and given by: (9) where is the intensity of photons reflected by both the object and the scattering media, and . is the intensity of background photons reflected by the scattering media. Then, the intensity of signal photons can be calculated by − . . We capture 9 measurements with objects embedded in the polyethylene foam at different depths, ranging from 2cm to 10cm with an increment of 1cm, to obtain the intensity of all the photons. Corresponding to the different object depths, the total thickness of the foam varies from 9cm to 17cm. The measured temporal responses are shown as the orange lines in Supplementary Fig. 10. Then we take the object out and measure the temporal responses to obtain the Supplementary Fig. 9 The relationship between the SBR and the object depth in the polyethylene foam. The SBR of the measurements at different depths can be obtained using Supplementary Equation (9). The results indicate that the signal photons decrease with the increase of the object depth. Using our setup, the maximum imaging depth of our algorithm to reconstruct the Lambertian object is 8cm where the SBR of the measurement is only 0.75%. background intensity, the blue lines in Supplementary Fig. 10. It can be seen that the difference between the two curves decreases with the increase in the object depth, which indicates the decrease of the signal photons. Using Supplementary Equation (9), we calculate the SBR of each measurement and find that the SBR is a decreasing function of the object depth, as shown in Supplementary Fig. 9.
The experiments are conducted with Lambertian objects at different transport mean free paths (TMFPs) to evaluate the validity range of our algorithm (The scattering lengths described below are all calculated in one direction). As a comparison, the reconstructions of the time gating algorithm are also presented in Supplementary Table 2. Since the time gating depends on isolating the signal photons, it can only work at less than 9.5 TMFPs and fails to handle the cases where the signal photons are too weak to be separated, such as the cases at a depth of 4cm to 10cm shown Supplementary Fig. 10 Comparison of the temporal response captured in the polyethylene foam with and without objects. The measurements are captured with different object depths in the polyethylene foam, ranging from 2cm to 10cm with an increment of 1cm. The orange lines show the photons reflected by both object and scattering media, while the blue lines show the photons reflected only by the scattering media. The difference between them indicates the intensity of the signal photons reflected by the object. in Supplementary Fig. 10. In contrast, benefitting from the accurate modeling and efficient use of all the optical photons, our algorithm performs well even when the object is at a depth of 25.4 TMFPs where the SBR is only 0.75%, as marked by the red dot in Supplementary Fig. 9. More importantly, it is found that the upper bound of the validity range of our algorithm can be improved further when we improve the SBR by increasing the incident power of the laser, using more sensitive photon detectors or enhancing the albedo of the object surface. For example, we use retroreflective objects to increase the signal intensity and test the algorithm's limit again under this condition, as shown in Supplementary We compare the results of time gating and our algorithm for reconstructing a Lambertian letter 'T' at different depths, ranging from 2cm to 9cm with an increment of 1cm. It can be found that time gating can only reconstruct the object at a depth less than 9.5 TMFPs, while our algorithm works even at 25.4 TMFPs. signal intensity, time gating can reconstruct the object at a depth less than 12.7 TMFPs, while our algorithm successfully reconstructs the object at a depth of 28.5 TMFPs.
To summarize, we demonstrate that all our experiments are within the validity range of the diffusion equation. Additionally, the validity range of our algorithm is theoretically unbounded but related to the SBR of the measurements in practical experiments. Using our present experimental setup, the maximum depth of the proposed algorithm to reconstruct the Lambertian and retroreflective object is 25.4 TMFPs and 28.5 TMFPs, respectively, i.e., the round-trip optical paths are 50.8 TMFPs and 57.0 TMFPs, respectively.  Table 3. Comparison of time gating and our algorithm for reconstructing retroreflective objects. We compare the results of time gating and our algorithm for reconstructing a retroreflective object at different depths, ranging from 2cm to 10cm with an increment of 1cm. It can be found that time gating can only reconstruct the object at a depth less than 12.7 TMFPs, while our algorithm works even at 28.5 TMFPs.

Supplementary Note 4: Evaluation of the resolution in the polyethylene foam scattering scenario
In this section, we derive and evaluate the resolution of our system in the polyethylene foam scattering scenario. The resolution is the minimal distinguishable distance Δ between two points & and ' , as shown in Supplementary  Fig. 11a. When the scanning interval is small enough that all diffused spots can cover the entire scanning area, the key factor affecting the resolution is the diffusion effect caused by the scattering media. In our analysis, we neglect the difference in axial propagation between the two points, and thus the scattering media primarily affect the lateral diffusion of the light. The lateral resolution of the system is derived using the full width at half maximum (FWHM) criterion. Supplementary Fig. 11. Illustration of the resolution analysis. a For a given depth , the minimal distinguishable distance Δ between two points & and ' is determined by the diffusion effect. b A full width at half maximum of a single diffusion spot measured at a depth of 5cm. The subfigure is the spatial distribution of a single diffusion spot. c Two points can be resolved if their distance is larger than Δ . The subfigure is the spatial distribution of double diffusion spots. d Resolution is calculated by Supplementary Equation (11) and is a function of depth. The theoretical resolutions at the depths of 3cm, 4cm, 5cm are 3.3cm, 4.2cm and 5.0cm, respectively. For a given depth , the propagation of the scattered light ( , , ) can be calculated by the theoretical solution of diffusion equation as shown in Supplementary Equation (1). Then the spatial distribution of the laterally diffused spot ^( , ) is given by: (10) where is the length of the time bins. The subfigure in Supplementary Fig. 11b shows the spatial distribution of diffused spots at a depth of 5cm, and the blue curve is the intensity along the central axis. Then the FWHM is defined as the width at the half maximum of the curve. As shown in Supplementary Fig. 11c, two points can be resolved from the measurement if the distance between them is larger than the FWHM. Therefore, the lateral resolution Δ of our system is given by: (11) where is the number of vertical pixels of the diffused spot shown in the subfigures of Supplementary Fig. 11b-c. Supplementary Fig. 11d illustrates that the resolution of the system calculated by Supplementary Equation (11) is an increasing function of the depth.
To verify the derived theoretical resolution of our system, we perform experiments based on the measurements of a resolution chart. There are four groups of stripe patterns with centers separated by 5cm, 4cm, 3cm, and 2cm, as shown in Supplementary Fig. 12a. Supplementary Fig. 12b-d shows the reconstructions of the chart which is embedded in the polyethylene foam at the depths of 3cm, 4cm, and 5cm, respectively. Correspondingly, the total thicknesses of the foam are 10cm, 11cm, and 12cm. The laser scans a grid of 64 × 64 points on the front surface of the media over a 45cm × 45cm square area so that the scanning interval is 0.7cm which is smaller than the theoretical resolution at all depths used in our experiments. According to Supplementary Fig. 11d, the theoretical resolutions at the depths of 3cm, 4cm, 5cm are 3.3cm, 4.2cm, and 5.0cm, respectively. As can be seen from Supplementary Fig.  12b, the patterns separated by 3cm are just exceed the predicted resolution at a depth of 3cm, while the groups with 4cm and 5cm intervals are clearly resolvable. For the result at a depth of 4cm shown in Supplementary Fig. 12c, the groups separated by 5cm are clearly discernible, and the groups with 4cm intervals are also resolvable but a little fuzzy because this interval is just at the boundary of the theoretical resolution. Supplementary Fig. 12d shows the result at a depth of 5cm, the patterns separated by 5cm are resolvable as predicted by the theoretical values. The comparisons of theoretical and measured resolutions at different depths are summarized in the table below Supplementary Fig. 12, which presents high consistency between the experimental results and the derived theoretical resolutions. In this section, we show the performance of our algorithm to reconstruct objects with different albedo surfaces. The difference in properties of the target and the scattering media directly leads to the difference of the photons returned by the targets and the scattering media, which further affects the reconstruction performance. We have demonstrated the ability of our algorithm in reconstructing Lambertian objects in the experiments above. It is extremely challenging since the Lambertian surfaces diffuse the photons isotropically so that the photons returned by the object are highly coupled with those returned by the scattering media, as shown in Supplementary Fig. 10. While, the retroreflective surfaces on the non-Lambertian objects can considerably increase the direct-reflected signal photons with a much larger difference from the background scattered ones, which aids in the reconstruction. Therefore, many works reconstruct the retroreflective objects. Here, we also evaluate the performance of our algorithm on such objects to demonstrate its robustness. Two-depth letters 'LU' and 'ST' are tested in the polyethylene foam scattering scenario in this experiment. As shown in Supplementary Table 4, the same letters are compared in three groups: both letters are Lambertian; one letter is Lambertian and the other is retroreflective; and the two letters are retroreflective. The depths of 'L' and 'U' are 3cm and 5cm, and the total thickness of the foam is 12cm. The depth of 'S' and 'T' are 4cm and 6cm, and the total thickness of the foam is 13cm. As a comparison, the results of the time gating algorithm are also given in Supplementary Table 4. It can be found that retroreflective surface allows for the better time-gated results of the closer objects such as 'U', while the time gating still does not work for the farther objects. In contrast, the proposed algorithm provides legible results in all three cases, which shows the robustness of our algorithm for reconstructing objects with different albedo surfaces. Additionally, we have shown that our algorithm can reconstruct the retroreflective object at the one-way optical depth of 28.5 TMFPs in Supplementary Table 3.

Reference
Captured Data Time gating Ours

Supplementary Table 4. Reconstructions of objects with different albedo surfaces.
Each set of two-depth letters is compared in three groups: both letters are Lambertian; one letter is Lambertian and the other is retroreflective; and the two letters are retroreflective.

Supplementary Note 6: Reconstructions of objects under ambient illumination
Different from the experiments above that are conducted in the dark environment, we demonstrate the ability of our algorithm under ambient illumination in this section. As shown in Supplementary Fig. 13a, 72W fluorescent lamps hanging 3m above the imaging scenario provide strong environmental lighting. Ambient illumination directly results in the increase of the background photons as the red curve shown in Supplementary Fig.13b, which makes the reconstruction more challenging. We compare the reconstructions of the two-depth letters 'ST' with different albedo surfaces under the ambient illumination, as shown in Supplementary Table 5. The depths of the 'S' and 'T' are 4cm and 6cm, respectively. The total thickness of the polyethylene foam is 13cm.
The change of illumination mode corresponds to the change of the source term on the right side of Eq. (1) in the main text. In our derivation, the source term is implicitly contained in two boundary conditions of the scattered field, Φ # " $ , % , & % and Φ " ( ! , " , ) in Eq. (4) and Eq. (6) in the main text. Since we use a boundary migration model of the scattered field to tackle the scattering problem rather than directly calculate the analytical solution of Eq. (1), the specific form of the source term is not important. As can be seen from Supplementary Table 5 In this section, we compare our algorithm with the methods f-k 14 and CDT 15 which are also based on the idea of boundary value migration and are proposed in seeing around the corner and seeing through the scattering media, respectively. Both imaging scenarios are shown in Supplementary Fig. 14a-b. The difference among these three methods in principle is that f-k is based on the wave equation which only considers the propagation of light in free space, CDT models the scattering effect as a blur kernel acting on the free-space propagation, while our algorithm accurately models the light in scattering media based on the diffusion equation.
The difference is also reflected in the temporal response, as shown in Supplementary Fig. 14. It can be found that the temporal responses captured in the NLoS scenario are the same as those in the ideal single scattering scenario as shown in Fig. 1a in the main text and carry the shape information of the object in the discontinuities of the transient measurements. According to Fermat's principle 15 , the types of transient discontinuities in temporal response, including local minimum, local maximum, or saddle point, correspond to the object's surface shape: convex, concave, or saddle, respectively. In this case, the red and green discontinuities of the temporal response in Supplementary Fig.  14 indicate that the vase surface shape is composed of a convex body and a saddle neck. With the increase in the scattering strength and the scattering region, the signal photons that carry the target information are drastically attenuated, as shown in Supplemental Fig. 14b or even almost disappeared in the volumetric scattering scenarios, as shown in Supplemental Fig. 14c. The temporal responses indicate that f-k and CDT both rely on the signal photons retained by light propagation in free space.
The reconstruction performance of the time gating, f-k, CDT and our algorithm are compared in the volumetric scattering scenario. First, we test the single Lambertian letter 'T' at different depths in the polyethylene foam, as shown in Supplementary Table 6. The temporal responses measured in such cases (Fig. 3b in the main text) show  that it is difficult to directly decouple the signal photons from the measurements. Therefore, time gating, f-k and CDT cannot recover the object information from non-gated measurements. After cropping the measurement according to the object depth prior, time gating recovers the object at a depth of 2cm, f-k and CDT recover the object at 2cm and 4cm. Since the CDT models and inverts the scattering process, it performs better than f-k. However, time gating, f-k and CDT are all relied on the signal photons and thus they cannot handle the case where the signal photons are dramatically attenuated by the scattering media. In contrast, our algorithm models the propagation of light in the scattering media and computationally uses all the photons, which allows us to recover the object at the depth from 2cm to 8cm using non-gated measurements and provides the best performance.
Then we compare the reconstructions of multi-depth 2D objects and 3D objects embedded in the polyethylene foam in Supplementary Table 7. Since these objects have complex depth distributions, signal photons reflected from different depths will overlap in their temporal distribution. In this case, even using the measurements taken by combining multiple time-shift gates corresponding to different depths of the object, time gating, f-k and CDT provide almost no discernible information to identify the signatures of the object. Although these three algorithms reconstruct the letter 'U' at a relatively close depth of 3cm, the signal-to-noise ratio (SNR) of the results is low. In contrast, our method correctly recovers all the targets with high SNR and also provides the best qualitative agreement with the ground truth.
Supplementary Table 8 shows the comparison of these algorithms in reconstructing 3D objects with complex structure in the fog chamber. Since the fog is a typical forward scattering media where backscattered light has little influence on signal photons, the photons reflected by the object and those do not reach the object can be distinguished in the temporal response taking depth as a priori, as shown in Fig. 4b in the main text. In this case, time gating and fk also provide partial discernible features of objects using cropped measurements. However, affected by the time variability of the fog, the temporal distribution of the signal photons is superimposed with dynamic noise so that time gating and f-k cannot reconstruct the exact shape of the object. Additionally, the CDT algorithm, which simply regards the scattering effect as a fuzzy kernel, is no longer suitable for the dynamic volumetric scattering media and thus fails to recover the objects embedded in the fog. In contrast, our algorithm successfully recovers all the objects, especially the details of objects such as the handle and lid of the 'Teapot' and the ear of the 'Bunny'.
Exposure time per scan position, photon counts per histogram and runtime of the above experiments are summarized in Supplementary Table 9. The runtime is measured on a MacBook Pro (2.3GHz Intel Core i5) under MATLAB environment. Even though we do not crop the data to remove the photons directly returned by the scattering media, the number of photons per histogram used in our experiment is 1-2 orders smaller than that in the CDT experiment as described in the respective paper. The average exposure time per scan position in our experiment is around 15 times shorter than that in the CDT experiment. Additionally, the runtime of our algorithm is almost equal to that of f-k and around 3 times faster than that of CDT. The runtime and the requirement of data acquisition including photon counts and exposure time show the potential of our algorithm for real-time applications.

Supplementary Note 8: Discussion of Depth Reconstruction
In this section, we discuss the ability of our algorithm in depth reconstruction. During our derivation, the depth information is included in the temporal boundary condition of the scattered field ( , , , = 0), and transformed into the frequency domain as Φ′( / , 0 , 1 ). To make it computationally feasible to establish the transformation relationship between the temporal and spatial boundary condition of the scattered field. The dispersion relation is used to transform the wave vector of depth 1 to frequency , so that Φ′( / , 0 , 1 ) is weighted and interpolated to be Φ c ( / , 0 , ). Then the spatial boundary condition ( , , = 0, ) is obtained by applying the inverse Fourier transform over / , 0 and a numerical transformation over using Eq. (6) in the main text. In this case, the depth information of the object is reflected in the temporal information of the measurements. However, the transformation from ( , , = 0, ) to Φ c ( / , 0 , ) over is an ill-conditioned problem as shown in Eq. (8) in Methods.
Specifically, the transfer matrix between and is a sparse matrix. Inverting sparse matrices is an ill-conditioned problem. Methods such as the Quasi-minimal residual method, Preconditioned Conjugate Gradient (PCG), or Generalized Minimal Residual (GMRES) can reduce the residuals of sparse matrix inversion at some level, while solving the inverse problem still affects the accuracy of our algorithm in depth reconstruction. We compare the depth reconstructions of time gating, f-k, CDT and our algorithm in Supplementary Table 10. The target is a bunny embedded in the fog chamber and the experimental condition is the same as that described in the main text. The front view, top view, and side view are given by the maximum projection of the 3D reconstructions in x-y, x-z, and y-z dimensions, respectively. The reconstructed depth information can be seen from the top and side view. The results captured without scattering media are used as the ground truth in this comparison. Since the timegated result is obtained by directly cropping the time bins according to the object depth prior, it reflects the depth distribution of the object to some extent. Suffering from the facts that f-k does not consider the scattering effect, CDT does not accurately model the propagation in the volumetric scattering media, and our algorithm is affected by the error in solving the ill-conditioned problem, these three methods can not accurately recover the depth of the object. Fortunately, the rough depth of the object can be directly estimated from the temporal response in the fog since the backscattered light has little effect on the objects, as shown in Fig. 4b in the main text. This estimated depth range can be used as a criterion to calibrate the original depth reconstructed by f-k, CDT and our algorithm. In this case, the top view and side view of the final reconstructions using these three methods can give the rough depth of the object as shown in the bottom three rows in Supplementary Table 10. Additionally, the depth reconstruction performance of fk is better than that of CDT and our algorithm because it recovers not only the rough depth range but also the details of the depth distribution of the object. Our algorithm loses part of depth information in the process of inverting the illconditioned problem and thus fails to reconstruct the accurate details of the depth distribution of objects. We will explore other relative algorithms aiming at solving the inversion of the large sparse matrix to help accurately recover the depth information in the future. The result captured without scattering media is used as the ground truth. The front view, top view, and side view are given by the maximum projection of the 3D reconstructions in x-y, x-z, and y-z dimensions, respectively.